import math
import numpy as np
import matplotlib.pyplot as plt
from itertools import permutations, combinations
from collections import Counter
import scipy.stats as stats


# 设置中文字体
def set_chinese_font():
    """设置中文字体"""
    plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
    plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

set_chinese_font()

def combinatorics_and_probability():
    print("=== 排列组合与概率统计 ===")
    
    # 1. 排列组合基础
    print("\n1. 排列组合基础")
    items = ['A', 'B', 'C', 'D']
    
    # 排列
    print("排列 P(4, 2):")
    perms = list(permutations(items, 2))
    for i, perm in enumerate(perms[:8], 1):
        print(f"  {i}. {perm}")
    print(f"  总排列数: {len(perms)}")
    print(f"  公式验证: P(4,2) = 4!/(4-2)! = {math.perm(4, 2)}")
    
    # 组合
    print("\n组合 C(4, 2):")
    combs = list(combinations(items, 2))
    for i, comb in enumerate(combs):
        print(f"  {i+1}. {comb}")
    print(f"  总组合数: {len(combs)}")
    print(f"  公式验证: C(4,2) = 4!/(2!×(4-2)!) = {math.comb(4, 2)}")
    
    # 2. 概率计算
    print("\n2. 概率计算")
    
    def binomial_probability(n, k, p):
        """二项分布概率"""
        comb = math.comb(n, k)
        return comb * (p ** k) * ((1 - p) ** (n - k))
    
    # 掷硬币概率
    n_coin = 10
    p_head = 0.5
    
    print(f"掷硬币{n_coin}次概率分布:")
    probabilities = []
    for k in range(n_coin + 1):
        prob = binomial_probability(n_coin, k, p_head)
        probabilities.append(prob)
        print(f"  正面出现{k}次: {prob:.4f} ({prob*100:.2f}%)")
    
    # 3. 生日悖论
    print("\n3. 生日悖论")
    
    def birthday_probability(n):
        """计算n个人中至少两人生日相同的概率"""
        if n > 365:
            return 1.0
        prob_no_match = 1.0
        for i in range(n):
            prob_no_match *= (365 - i) / 365
        return 1 - prob_no_match
    
    for n in [23, 30, 40, 50]:
        prob = birthday_probability(n)
        print(f"  {n}人中至少两人生日相同的概率: {prob:.4f} ({prob*100:.2f}%)")
    
    # 4. 统计描述
    print("\n4. 统计描述分析")
    
    # 生成模拟数据（身高数据）
    np.random.seed(42)  # 设置随机种子保证结果可重现
    data = np.random.normal(170, 10, 1000)  # 均值170，标准差10
    
    print(f"数据统计描述 (1000个样本):")
    print(f"  均值: {np.mean(data):.2f}")
    print(f"  中位数: {np.median(data):.2f}")
    
    # 修复众数计算
    try:
        # 方法1: 使用scipy.stats.mode
        mode_result = stats.mode(data)
        # 不同版本的scipy返回格式不同，需要兼容处理
        if hasattr(mode_result, 'mode'):
            mode_value = mode_result.mode[0] if hasattr(mode_result.mode, '__getitem__') else mode_result.mode
        else:
            mode_value = mode_result[0]
        print(f"  众数: {mode_value:.2f}")
    except:
        # 方法2: 手动计算众数
        hist, bin_edges = np.histogram(data, bins=30)
        max_bin_index = np.argmax(hist)
        mode_value = (bin_edges[max_bin_index] + bin_edges[max_bin_index + 1]) / 2
        print(f"  众数(近似): {mode_value:.2f}")
    
    print(f"  标准差: {np.std(data):.2f}")
    print(f"  方差: {np.var(data):.2f}")
    print(f"  极差: {np.ptp(data):.2f}")
    print(f"  四分位距: {np.percentile(data, 75) - np.percentile(data, 25):.2f}")
    
    # 5. 条件概率和贝叶斯定理
    print("\n5. 条件概率与贝叶斯定理")
    
    # 疾病检测例子
    prevalence = 0.01  # 疾病患病率 1%
    sensitivity = 0.95  # 检测灵敏度 95%
    specificity = 0.90  # 检测特异度 90%
    
    # 计算检测阳性的情况下真正患病的概率
    p_disease = prevalence
    p_positive_given_disease = sensitivity
    p_positive_given_no_disease = 1 - specificity
    
    p_positive = (p_disease * p_positive_given_disease + 
                 (1 - p_disease) * p_positive_given_no_disease)
    
    p_disease_given_positive = (p_disease * p_positive_given_disease) / p_positive
    
    print(f"疾病检测贝叶斯分析:")
    print(f"  患病率: {p_disease:.3f}")
    print(f"  检测灵敏度: {sensitivity:.3f}")
    print(f"  检测特异度: {specificity:.3f}")
    print(f"  检测阳性时真正患病的概率: {p_disease_given_positive:.3f}")
    
    # 可视化
    visualize_probability_stats(probabilities, data, n_coin)

def visualize_probability_stats(probabilities, data, n_coin):
    """可视化概率统计结果"""
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    
    # 1. 二项分布
    k_values = range(n_coin + 1)
    axes[0, 0].bar(k_values, probabilities, alpha=0.7, color='skyblue', edgecolor='black')
    axes[0, 0].set_xlabel('成功次数 (k)')
    axes[0, 0].set_ylabel('概率')
    axes[0, 0].set_title('二项分布 B(10, 0.5)')
    axes[0, 0].grid(True, alpha=0.3)
    
    # 添加概率值标签
    for k, prob in zip(k_values, probabilities):
        if prob > 0.05:  # 只显示概率较大的标签
            axes[0, 0].text(k, prob + 0.01, f'{prob:.3f}', 
                           ha='center', va='bottom', fontsize=8)
    
    # 2. 生日悖论
    n_people = range(1, 61)
    birthday_probs = [birthday_probability(n) for n in n_people]
    axes[0, 1].plot(n_people, birthday_probs, 'r-', linewidth=2)
    axes[0, 1].axhline(y=0.5, color='gray', linestyle='--', alpha=0.7)
    axes[0, 1].axvline(x=23, color='gray', linestyle='--', alpha=0.7)
    axes[0, 1].set_xlabel('人数')
    axes[0, 1].set_ylabel('概率')
    axes[0, 1].set_title('生日悖论')
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].text(23, 0.5, '23人 → 50%', ha='right', va='bottom')
    
    # 3. 数据直方图
    axes[0, 2].hist(data, bins=30, density=True, alpha=0.7, color='lightgreen', edgecolor='black')
    axes[0, 2].set_xlabel('身高 (cm)')
    axes[0, 2].set_ylabel('密度')
    axes[0, 2].set_title('身高数据分布')
    axes[0, 2].grid(True, alpha=0.3)
    
    # 添加正态分布曲线
    x = np.linspace(data.min(), data.max(), 100)
    pdf = stats.norm.pdf(x, data.mean(), data.std())
    axes[0, 2].plot(x, pdf, 'r-', linewidth=2, label='正态分布拟合')
    axes[0, 2].legend()
    
    # 4. 箱线图
    axes[1, 0].boxplot(data, vert=True, patch_artist=True)
    axes[1, 0].set_title('数据箱线图')
    axes[1, 0].set_ylabel('身高 (cm)')
    axes[1, 0].grid(True, alpha=0.3)
    
    # 5. 累积分布函数
    data_sorted = np.sort(data)
    cdf = np.arange(1, len(data) + 1) / len(data)
    axes[1, 1].plot(data_sorted, cdf, 'b-', linewidth=2)
    axes[1, 1].set_xlabel('身高 (cm)')
    axes[1, 1].set_ylabel('累积概率')
    axes[1, 1].set_title('累积分布函数 (CDF)')
    axes[1, 1].grid(True, alpha=0.3)
    
    # 6. QQ图（正态性检验）
    stats.probplot(data, dist="norm", plot=axes[1, 2])
    axes[1, 2].set_title('QQ图 - 正态性检验')
    
    plt.tight_layout()
    plt.show()

def advanced_probability_examples():
    """高级概率例子"""
    print("\n=== 高级概率例子 ===")
    
    # 1. 蒙提霍尔问题
    print("\n1. 蒙提霍尔问题 (三门问题)")
    
    def monty_hall_simulation(n_trials=10000):
        stay_wins = 0
        switch_wins = 0
        
        for _ in range(n_trials):
            # 随机放置汽车
            doors = [0, 0, 0]  # 0代表山羊，1代表汽车
            car_position = np.random.randint(0, 3)
            doors[car_position] = 1
            
            # 玩家随机选择一扇门
            player_choice = np.random.randint(0, 3)
            
            # 主持人打开一扇有山羊的门（不是玩家选择的，也不是有汽车的）
            host_choices = [i for i in range(3) if i != player_choice and doors[i] == 0]
            host_choice = np.random.choice(host_choices)
            
            # 计算坚持和换门的获胜次数
            if doors[player_choice] == 1:
                stay_wins += 1
            
            # 换到另一扇未打开的门
            switch_choice = [i for i in range(3) if i != player_choice and i != host_choice][0]
            if doors[switch_choice] == 1:
                switch_wins += 1
        
        return stay_wins / n_trials, switch_wins / n_trials
    
    stay_prob, switch_prob = monty_hall_simulation(10000)
    print(f"  坚持原选择的获胜概率: {stay_prob:.4f}")
    print(f"  换门的获胜概率: {switch_prob:.4f}")
    print(f"  理论值: 坚持=1/3≈0.333, 换门=2/3≈0.667")
    
    # 2. 泊松分布
    print("\n2. 泊松分布例子")
    
    def poisson_probability(lambd, k):
        """泊松分布概率"""
        return (lambd ** k * math.exp(-lambd)) / math.factorial(k)
    
    # 模拟商店每小时顾客数
    lambd = 5  # 平均每小时5个顾客
    print(f"平均每小时{lambd}个顾客的概率分布:")
    for k in range(10):
        prob = poisson_probability(lambd, k)
        print(f"  {k}个顾客: {prob:.4f}")
    
    # 3. 几何分布
    print("\n3. 几何分布例子")
    
    def geometric_probability(p, k):
        """几何分布概率"""
        return (1 - p) ** (k - 1) * p
    
    # 抛硬币直到出现正面
    p_head = 0.5
    print(f"抛硬币直到出现正面的概率分布:")
    for k in range(1, 6):
        prob = geometric_probability(p_head, k)
        print(f"  第{k}次出现正面: {prob:.4f}")

def birthday_probability(n):
    """计算n个人中至少两人生日相同的概率"""
    if n > 365:
        return 1.0
    prob_no_match = 1.0
    for i in range(n):
        prob_no_match *= (365 - i) / 365
    return 1 - prob_no_match

def practical_applications():
    """实际应用案例"""
    print("\n=== 实际应用案例 ===")
    
    # 1. 彩票中奖概率
    print("\n1. 双色球中奖概率计算")
    
    def lottery_probability():
        # 双色球规则: 33选6 + 16选1
        total_red = math.comb(33, 6)
        total_blue = 16
        total_combinations = total_red * total_blue
        
        print(f"总组合数: {total_combinations:,}")
        print("各奖项中奖概率:")
        
        # 一等奖: 6+1
        prob_1st = 1 / total_combinations
        print(f"  一等奖: 1/{total_combinations:,} ≈ {prob_1st:.10f}")
        
        # 二等奖: 6+0
        prob_2nd = (math.comb(6, 6) * math.comb(27, 0) * 15) / total_combinations
        print(f"  二等奖: {prob_2nd:.8f}")
        
        # 三等奖: 5+1
        prob_3rd = (math.comb(6, 5) * math.comb(27, 1) * 1) / total_combinations
        print(f"  三等奖: {prob_3rd:.8f}")
    
    lottery_probability()
    
    # 2. 产品质量检验
    print("\n2. 产品质量检验抽样方案")
    
    def quality_control():
        n_samples = 50  # 抽样数量
        defect_rate = 0.02  # 缺陷率
        acceptance_number = 1  # 可接受缺陷数
        
        # 计算接受概率
        prob_accept = 0
        for k in range(acceptance_number + 1):
            prob = binomial_probability(n_samples, k, defect_rate)
            prob_accept += prob
        
        print(f"抽样方案: n={n_samples}, c={acceptance_number}")
        print(f"缺陷率: {defect_rate*100}%")
        print(f"接受概率: {prob_accept:.4f}")
        
        # 绘制OC曲线
        defect_rates = np.linspace(0, 0.1, 100)
        accept_probs = []
        
        for p in defect_rates:
            prob_acc = 0
            for k in range(acceptance_number + 1):
                prob_acc += binomial_probability(n_samples, k, p)
            accept_probs.append(prob_acc)
        
        plt.figure(figsize=(10, 6))
        plt.plot(defect_rates, accept_probs, 'b-', linewidth=2)
        plt.axvline(x=defect_rate, color='r', linestyle='--', label=f'当前缺陷率 ({defect_rate*100}%)')
        plt.xlabel('缺陷率')
        plt.ylabel('接受概率')
        plt.title('操作特性曲线 (OC曲线)')
        plt.grid(True, alpha=0.3)
        plt.legend()
        plt.show()
    
    quality_control()

def binomial_probability(n, k, p):
    """二项分布概率"""
    comb = math.comb(n, k)
    return comb * (p ** k) * ((1 - p) ** (n - k))

# 运行所有函数
if __name__ == "__main__":
    combinatorics_and_probability()
    advanced_probability_examples()
    practical_applications()